//-------------------------------------------------------------------------------------------------
//  matrix4x4.h - 4x4 square matrix template implementation
//  (c) 2009 Simon Miessler
//-------------------------------------------------------------------------------------------------

#pragma once

#include "math_exports.hpp"
#include "vector4.hpp"
#include "matrix3x3.hpp"

#pragma pack(push)
#pragma pack(1)




namespace math
{
	/**
	 * @brief This class represents a 4x4 square, row-major matrix.
	 *
	 * This matrix can be seen as a rotation and translation matrix in R3.
	 * It's the only matrix from the math library, that has a DirectX counterpart: D3DXMATRIX.
	 * Every row represents an axis of the new coordinate system:
	 * - first row:  x-axis
	 * - second row: y-axis
	 * - third row:  z-axis
	 * - fourth row: translation vector
	 */
	template
	<
		typename T,
		typename eq_pred = equality_predicate<T, equality_digits<T>::digits>
	>
	class matrix4x4
		: boost::equality_comparable< matrix4x4<T, eq_pred>
		, boost::addable< matrix4x4<T, eq_pred>
		, boost::subtractable< matrix4x4<T, eq_pred>
		, boost::dividable2< matrix4x4<T, eq_pred>, T
		, boost::multipliable< matrix4x4<T, eq_pred>
		, boost::multipliable2< matrix4x4<T, eq_pred>, T
		> > > > > >
	{
	public:

		/// The used floating-point type
		typedef T   float_t;

	public:

		/// First column, first row
		T _11;
		/// Second column, first row
		T _12;
		/// Third column, first row
		T _13;
		/// Fourth column, first row
		T _14;
		/// First column, second row
		T _21;
		/// Second column, second row
		T _22;
		/// Third column, second row
		T _23;
		/// Fourth column, second row
		T _24;
		/// First column, third row
		T _31;
		/// Second column, third row
		T _32;
		/// Third column, third row
		T _33;
		/// Fourth column, third row
		T _34;
		/// First column, fourth row
		T _41;
		/// Second column, fourth row
		T _42;
		/// Third column, fourth row
		T _43;
		/// Fourth column, fourth row
		T _44;


		/**
		 * Create a matrix. The matrix will be an identity matrix by default.
		 */
		matrix4x4()
		{
			identity();
		}

		/**
		 * Create a matrix from 4 vectors.
		 * This function can be used when this matrix is used as a translation
		 * matrix in R3.
		 *
		 * @param xAxis      the new coordinate system's x-axis
		 * @param yAxis      the new coordinate system's y-axis
		 * @param zAxis      the new coordinate system's z-axis
		 * @param position   the new coordinate system's position
		 */
		matrix4x4(const vector3<T, eq_pred>& xAxis,
				  const vector3<T, eq_pred>& yAxis,
				  const vector3<T, eq_pred>& zAxis,
				  const vector3<T, eq_pred>& position)
			: _11(xAxis.x),    _12(xAxis.y),    _13(xAxis.z),    _14(0)
			, _21(yAxis.x),    _22(yAxis.y),    _23(yAxis.z),    _24(0)
			, _31(zAxis.x),    _32(zAxis.y),    _33(zAxis.z),    _34(0)
			, _41(position.x), _42(position.y), _43(position.z), _44(1)
		{}

		/**
		 * Create a matrix from 4 vectors.
		 */
		matrix4x4(const vector4<T, eq_pred>& column1,
				  const vector4<T, eq_pred>& column2,
				  const vector4<T, eq_pred>& column3,
				  const vector4<T, eq_pred>& column4)
			: _11(column1.x), _12(column1.y), _13(column1.z), _14(column1.w)
			, _21(column2.x), _22(column2.y), _23(column2.z), _24(column2.w)
			, _31(column3.x), _32(column3.y), _33(column3.z), _34(column3.w)
			, _41(column4.x), _42(column4.y), _43(column4.z), _44(column4.w)
		{}

		/**
		 * Create a matrix from a 3x3 matrix and a position vector.
		 * This function can be used when this matrix is used as a translation
		 * matrix in R3.
		 *
		 * @param axes       the new coordinate system's axes (x, y and z)
		 * @param position   the new coordinate system's position
		 */
		matrix4x4(const matrix3x3<T, eq_pred>& axes,
				  const vector3<T, eq_pred>& position)
			: _11(axes._11),   _12(axes._12),   _13(axes._13),   _14(0)
			, _21(axes._21),   _22(axes._22),   _23(axes._23),   _24(0)
			, _31(axes._31),   _32(axes._32),   _33(axes._33),   _34(0)
			, _41(position.x), _42(position.y), _43(position.z), _44(1)
		{}

		/**
		 * Create a matrix4x4 from scratch.
		 * Every four values form a column.
		 */
		matrix4x4(T _11, T _12, T _13, T _14,
				  T _21, T _22, T _23, T _24,
				  T _31, T _32, T _33, T _34,
				  T _41, T _42, T _43, T _44)
			: _11(_11), _12(_12), _13(_13), _14(_14)
			, _21(_21), _22(_22), _23(_23), _24(_24)
			, _31(_31), _32(_32), _33(_33), _34(_34)
			, _41(_41), _42(_42), _43(_43), _44(_44)
		{}

		/**
		 * Create a matrix4x4 from another matrix, using a different type.
		 */
		template
		<
			typename Y,
			typename eq_pred
		>
		explicit matrix4x4(const matrix4x4<Y, eq_pred>& matrix)
			: _11((T)matrix._11), _12((T)matrix._12), _13((T)matrix._13), _14((T)matrix._14)
			, _21((T)matrix._21), _22((T)matrix._22), _23((T)matrix._23), _24((T)matrix._24)
			, _31((T)matrix._31), _32((T)matrix._32), _33((T)matrix._33), _34((T)matrix._34)
			, _41((T)matrix._41), _42((T)matrix._42), _43((T)matrix._43), _44((T)matrix._44)
		{}



		/**
		 * Set this matrix to an identity matrix.
		 *
		 * @return a reference to this matrix
		 */
		matrix4x4<T, eq_pred>& identity()
		{
			_11 = 1; _12 = 0; _13 = 0; _14 = 0;
			_21 = 0; _22 = 1; _23 = 0; _24 = 0;
			_31 = 0; _32 = 0; _33 = 1; _34 = 0;
			_41 = 0; _42 = 0; _43 = 0; _44 = 1;

			return *this;
		}

		/**
		 * Test if this matrix is an identity matrix.
		 *
		 * @return true    when this matrix is an identity matrix
		 * @return false   when it is not
		 */
		bool isIdentity() const
		{
			// Instead of using operator== directly, we use the equality_predicate
			// that is specialized for certain types (floating-point for the most)
			return row1() == vector4<T, eq_pred>(1, 0, 0, 0) &&
				   row2() == vector4<T, eq_pred>(0, 1, 0, 0) &&
				   row3() == vector4<T, eq_pred>(0, 0, 1, 0) &&
				   row4() == vector4<T, eq_pred>(0, 0, 0, 1);
		}

		/**
		 * Calculate the determinant of this matrix.
		 *
		 * @see   http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
		 */
		T determinant() const
		{
			return _11 * _22 * _33 * _44 + _11 * _23 * _34 * _42 + _11 * _24 * _32 * _43
				 + _12 * _21 * _34 * _43 + _12 * _23 * _31 * _44 + _12 * _24 * _33 * _41
				 + _13 * _21 * _32 * _44 + _13 * _22 * _34 * _41 + _13 * _24 * _31 * _42
				 + _14 * _21 * _33 * _42 + _14 * _22 * _31 * _43 + _14 * _23 * _32 * _41
				 - _11 * _22 * _34 * _43 - _11 * _23 * _32 * _44 - _11 * _24 * _33 * _42
				 - _12 * _21 * _33 * _44 - _12 * _23 * _34 * _41 - _12 * _24 * _31 * _43
				 - _13 * _21 * _34 * _42 - _13 * _22 * _31 * _44 - _13 * _24 * _32 * _41
				 - _14 * _21 * _32 * _43 - _14 * _22 * _33 * _41 - _14 * _23 * _31 * _42;
		}

		/**
		 * Set this matrix to its inverse matrix.
		 *
		 * @return a reference to this matrix
		 * @see   http://www.cvl.iis.u-tokyo.ac.jp/~miyazaki/tech/teche23.html
		 */
		matrix4x4<T, eq_pred>& inverse()
		{
			T A = determinant();
			assert(A != 0.0f);

			*this = 1 / A * matrix4x4<T, eq_pred>(_22 * _33 * _44 + _23 * _34 * _42 + _24 * _32 * _43 - _22 * _34 * _43 - _23 * _32 * _44 - _24 * _33 * _42, // 11
												  _12 * _34 * _43 + _13 * _32 * _44 + _14 * _33 * _42 - _12 * _33 * _44 - _13 * _34 * _42 - _14 * _32 * _43, // 12
												  _12 * _23 * _44 + _13 * _24 * _42 + _14 * _22 * _43 - _12 * _24 * _43 - _13 * _22 * _44 - _14 * _23 * _42, // 13
												  _12 * _24 * _33 + _13 * _22 * _34 + _14 * _23 * _32 - _12 * _23 * _34 - _13 * _24 * _32 - _14 * _22 * _33, // 14
												  _21 * _34 * _43 + _23 * _31 * _44 + _24 * _33 * _41 - _21 * _33 * _44 - _23 * _34 * _41 - _24 * _31 * _43, // 21
												  _11 * _33 * _44 + _13 * _34 * _41 + _14 * _31 * _43 - _11 * _34 * _43 - _13 * _31 * _44 - _14 * _33 * _41, // 22
												  _11 * _24 * _43 + _13 * _21 * _44 + _14 * _23 * _41 - _11 * _23 * _44 - _13 * _24 * _41 - _14 * _21 * _43, // 23
												  _11 * _23 * _34 + _13 * _24 * _31 + _14 * _21 * _33 - _11 * _24 * _33 - _13 * _21 * _34 - _14 * _23 * _31, // 24
												  _21 * _32 * _44 + _22 * _34 * _41 + _24 * _31 * _42 - _21 * _34 * _42 - _22 * _31 * _44 - _24 * _32 * _41, // 31
												  _11 * _34 * _42 + _12 * _31 * _44 + _14 * _32 * _41 - _11 * _32 * _44 - _12 * _34 * _41 - _14 * _31 * _42, // 32
												  _11 * _22 * _44 + _12 * _24 * _41 + _14 * _21 * _42 - _11 * _24 * _42 - _12 * _21 * _44 - _14 * _22 * _41, // 33
												  _11 * _24 * _32 + _12 * _21 * _34 + _14 * _22 * _31 - _11 * _22 * _34 - _12 * _24 * _31 - _14 * _21 * _32, // 34
												  _21 * _33 * _42 + _22 * _31 * _43 + _23 * _32 * _41 - _21 * _32 * _43 - _22 * _33 * _41 - _23 * _31 * _42, // 41
												  _11 * _32 * _43 + _12 * _33 * _41 + _13 * _31 * _42 - _11 * _33 * _42 - _12 * _31 * _43 - _13 * _32 * _41, // 42
												  _11 * _23 * _42 + _12 * _21 * _43 + _13 * _22 * _41 - _11 * _22 * _43 - _12 * _23 * _41 - _13 * _21 * _42, // 43
												  _11 * _22 * _33 + _12 * _23 * _31 + _13 * _21 * _32 - _11 * _23 * _32 - _12 * _21 * _33 - _13 * _22 * _31  // 44
										);

			return *this;
		}

		/**
		 * Set this matrix to its transpose matrix.
		 *
		 * @return a reference to this matrix
		 */
		matrix4x4<T, eq_pred>& transpose()
		{
			std::swap(_12, _21);
			std::swap(_13, _31);
			std::swap(_14, _41);
			std::swap(_23, _32);
			std::swap(_24, _42);
			std::swap(_34, _43);

			return *this;
		}

		/**
		 * Get the rotation part of this matrix (represented by a 3x3 matrix).
		 */
		matrix3x3<T, eq_pred> rotation() const
		{
			matrix3x3<T, eq_pred> tmp(_11, _12, _13,
									  _21, _22, _23,
									  _31, _32, _33);

			return tmp;
		}



		/**
		 * Get the first row of the matrix.
		 */
		vector4<T, eq_pred> row1() const
		{
			vector4<T, eq_pred> tmp(_11, _12, _13, _14);

			return tmp;
		}

		/**
		 * Get the second row of the matrix.
		 */
		vector4<T, eq_pred> row2() const
		{
			vector4<T, eq_pred> tmp(_21, _22, _23, _24);

			return tmp;
		}

		/**
		 * Get the third row of the matrix.
		 */
		vector4<T, eq_pred> row3() const
		{
			vector4<T, eq_pred> tmp(_31, _32, _33, _34);

			return tmp;
		}

		/**
		 * Get the fourth row of the matrix.
		 */
		vector4<T, eq_pred> row4() const
		{
			vector4<T, eq_pred> tmp(_41, _42, _43, _44);

			return tmp;
		}



		/**
		 * Get the first column of the matrix.
		 */
		vector4<T, eq_pred> col1() const
		{
			vector4<T, eq_pred> tmp(_11, _21, _31, _41);

			return tmp;
		}

		/**
		 * Get the second column of the matrix.
		 */
		vector4<T, eq_pred> col2() const
		{
			vector4<T, eq_pred> tmp(_12, _22, _32, _42);

			return tmp;
		}

		/**
		 * Get the third column of the matrix.
		 */
		vector4<T, eq_pred> col3() const
		{
			vector4<T, eq_pred> tmp(_13, _23, _33, _43);

			return tmp;
		}

		/**
		 * Get the fourth column of the matrix.
		 */
		vector4<T, eq_pred> col4() const
		{
			vector4<T, eq_pred> tmp(_14, _24, _34, _44);

			return tmp;
		}



		/**
		 * Add another 4x4 matrix to this one.
		 */
		matrix4x4<T, eq_pred>& operator+=(const matrix4x4<T, eq_pred>& matrix)
		{
			_11 += matrix._11; _12 += matrix._12; _13 += matrix._13; _14 += matrix._14;
			_21 += matrix._21; _22 += matrix._22; _23 += matrix._23; _24 += matrix._24;
			_31 += matrix._31; _32 += matrix._32; _33 += matrix._33; _34 += matrix._34;
			_41 += matrix._41; _42 += matrix._42; _43 += matrix._43; _44 += matrix._44;

			return *this;
		}

		/**
		 * Subtract another 4x4 matrix from this one.
		 */
		matrix4x4<T, eq_pred>& operator-=(const matrix4x4<T, eq_pred>& mat)
		{
			_11 -= mat._11; _12 -= mat._12; _13 -= mat._13; _14 -= mat._14;
			_21 -= mat._21; _22 -= mat._22; _23 -= mat._23; _24 -= mat._24;
			_31 -= mat._31; _32 -= mat._32; _33 -= mat._33; _34 -= mat._34;
			_41 -= mat._41; _42 -= mat._42; _43 -= mat._43; _44 -= mat._44;

			return *this;
		}

		/**
		 * Multiply this 4x4 matrix by a scalar.
		 */
		matrix4x4<T, eq_pred>& operator*=(T scalar)
		{
			_11 *= scalar; _12 *= scalar; _13 *= scalar; _14 *= scalar;
			_21 *= scalar; _22 *= scalar; _23 *= scalar; _24 *= scalar;
			_31 *= scalar; _32 *= scalar; _33 *= scalar; _34 *= scalar;
			_41 *= scalar; _42 *= scalar; _43 *= scalar; _44 *= scalar;

			return *this;
		}

		/**
		 * Multiply this 4x4 matrix by another matrix.
		 */
		matrix4x4<T, eq_pred>& operator*=(const matrix4x4<T, eq_pred>& m)
		{
			//         Matrix                   Matrix        
			// ( 0,0  0,1  0,2  0,3 )   ( 0,0  0,1  0,2  0,3 )
			// ( 1,0  1,1  1,2  1,3 ) * ( 1,0  1,1  1,2  1,3 )
			// ( 2,0  2,1  2,2  2,3 )   ( 2,0  2,1  2,2  2,3 )
			// ( 3,0  3,1  3,2  3,3 )   ( 3,0  3,1  3,2  3,3 )
			*this = matrix4x4<T, eq_pred>(dot(row1(), m.col1()), dot(row1(), m.col2()), dot(row1(), m.col3()), dot(row1(), m.col4()),
										  dot(row2(), m.col1()), dot(row2(), m.col2()), dot(row2(), m.col3()), dot(row2(), m.col4()),
										  dot(row3(), m.col1()), dot(row3(), m.col2()), dot(row3(), m.col3()), dot(row3(), m.col4()),
										  dot(row4(), m.col1()), dot(row4(), m.col2()), dot(row4(), m.col3()), dot(row4(), m.col4()));

			return *this;
		}

		/**
		 * Divide this 4x4 matrix by a scalar.
		 */
		matrix4x4<T, eq_pred>& operator/=(T scalar)
		{
			_11 /= scalar; _12 /= scalar; _13 /= scalar; _14 /= scalar;
			_21 /= scalar; _22 /= scalar; _23 /= scalar; _24 /= scalar;
			_31 /= scalar; _32 /= scalar; _33 /= scalar; _34 /= scalar;
			_41 /= scalar; _42 /= scalar; _43 /= scalar; _44 /= scalar;

			return *this;
		}

		/**
		 * Test if this matrix is equal to another one.
		 */
		bool operator==(const matrix4x4<T, eq_pred>& that) const
		{
			// Instead of using operator== directly, we use the equality_predicate
			// that is specialized for certain types (floating-point for the most)
			eq_pred p;
			return p(_11, that._11) && p(_12, that._12) && p(_13, that._13) && p(_14, that._14) &&
				   p(_21, that._21) && p(_22, that._22) && p(_33, that._33) && p(_24, that._24) &&
				   p(_31, that._31) && p(_32, that._32) && p(_33, that._33) && p(_34, that._34) &&
				   p(_41, that._41) && p(_42, that._42) && p(_43, that._43) && p(_44, that._44);
		}



		/**
		 * Create a 4x4 scale matrix that scales in all three axes.
		 */
		static matrix4x4<T, eq_pred> scale(T scale)
		{
			matrix4x4<T, eq_pred> tmp(matrix3x3<T, eq_pred>::scale(scale),
									  vector3<T, eq_pred>(0, 0, 0));

			return tmp;
		}

		/**
		 * Create a 4x4 scale matrix that scales individually per axis.
		 */
		static matrix4x4<T, eq_pred> scale(T xScale, T yScale, T zScale)
		{
			matrix4x4<T, eq_pred> tmp(matrix3x3<T, eq_pred>::scale(xScale, yScale, zScale),
									  vector3<T, eq_pred>(0, 0, 0));

			return tmp;
		}

		/**
		 * Create a 4x4 scale matrix that scales individually per axis.
		 */
		static matrix4x4<T, eq_pred> scale(const vector3<T, eq_pred>& scale)
		{
			matrix4x4<T, eq_pred> tmp(matrix3x3<T, eq_pred>::scale(scale.x, scale.y, scale.z),
									  vector3<T, eq_pred>(0, 0, 0));

			return tmp;
		}



		/**
		 * Create a left-handed perspective projection matrix based on a field of view.
		 *
		 * @param fovY     the view-angle (horizontal) in radians
		 * @param aspect   the aspect ratio of the view
		 * @param zn       z-near distance from the camera in world units
		 * @param zf       z-far distance from the camera in world units
		 *
		 * @return   a left-handed perspective projection matrix
		 */
		static matrix4x4<T, eq_pred> perspectiveLH(angle<T> fovY, T aspect, T zn, T zf)
		{
			angle<T> fovYHalf = fovY / scalar<T>(2);

			T yScale = fovYHalf.cos() / fovYHalf.sin();
			T xScale = yScale / aspect;

			matrix4x4<T, eq_pred> tmp(xScale, 0,      0,              0,
									  0,      yScale, 0,              0,
									  0,      0,      zf/(zf-zn),     1,
									  0,      0,      -zn*zf/(zf-zn), 0);

			return tmp;
		}

		/**
		 * Create a right-handed perspective projection matrix based on a field of view.
		 *
		 * @param fovY     the view-angle (horizontal) in radians
		 * @param aspect   the aspect ratio of the view
		 * @param zn       z-near distance from the camera in world units
		 * @param zf       z-far distance from the camera in world units
		 *
		 * @return   a right-handed perspective projection matrix
		 */
		static matrix4x4<T, eq_pred> perspectiveRH(angle<T> fovY, T aspect, T zn, T zf)
		{
			angle<T> fovYHalf = fovY / scalar<T>(2);

			T yScale = fovYHalf.cos() / fovYHalf.sin();
			T xScale = yScale / aspect;

			matrix4x4<T, eq_pred> tmp(xScale, 0,      0,             0,
									  0,      yScale, 0,             0,
									  0,      0,      zf/(zn-zf),   -1,
									  0,      0,      zn*zf/(zn-zf), 0);

			return tmp;
		}



		/**
		 * Create a left-handed, orthogonal projection matrix.
		 *
		 * @param width    the width of the view volume
		 * @param height   the height of the view volume
		 * @param zn       z-near distance from the camera in world units
		 * @param zf       z-far distance from the camera in world units
		 *
		 * @return   a left-handed orthogonal projection matrix
		 */
		static matrix4x4<T, eq_pred> orthoLH(T width, T height, T zn, T zf)
		{
			matrix4x4<T, eq_pred> tmp(2/width, 0,        0,           0,
									  0,       2/height, 0,           0,
									  0,       0,        1/(zf-zn),   0,
									  0,       0,        -zn/(zf-zn), 1);

			return tmp;
		}

		/**
		 * Create a right-handed, orthogonal projection matrix.
		 *
		 * @param width    the width of the view volume
		 * @param height   the height of the view volume
		 * @param zn       z-near distance from the camera in world units
		 * @param zf       z-far distance from the camera in world units
		 *
		 * @return   a right-handed orthogonal projection matrix
		 */
		static matrix4x4<T, eq_pred> orthoRH(T width, T height, T zn, T zf)
		{
			matrix4x4<T, eq_pred> tmp(2/width, 0,        0,          0,
									  0,       2/height, 0,          0,
									  0,       0,        1/(zn-zf),  0,
									  0,       0,        zn/(zn-zf), 1);

			return tmp;
		}

		/**
		 * Create a left-handed, orthogonal off center projection matrix.
		 *
		 * @param left     the leftmost coordinate
		 * @param right    the rightmost coordinate
		 * @param top      the topmost coordinate
		 * @param bottom   the bottom most coordinate
		 * @param zn       z-near distance from the camera in world units
		 * @param zf       z-far distance from the camera in world units
		 */
		static matrix4x4<T, eq_pred> orthoOffCenterLH(T left, T right, T top, T bottom, T zNear, T zFar)
		{
			return matrix4x4<T, eq_pred>(2 / (right - left),              0,                               0,                      0,
										 0,                               2 / (top - bottom),              0,                      0,
										 0,                               0,                               1 / (zFar - zNear),     0,
										 (left + right) / (left - right), (top + bottom) / (bottom - top), zNear / (zNear - zFar), 1);
		}




		/**
		 * Create a left-handed, look-at matrix.
		 *
		 * @param pos   the position of the camera in world space
		 * @param at    the target of the camera in world space
		 * @param up    the up vector of the world, usually (0, 1, 0)
		 *
		 * @return   a left-handed, look-at matrix
		 */
		static matrix4x4<T, eq_pred> lookAtLH(const vector3<T, eq_pred>& pos, const vector3<T, eq_pred>& at, const vector3<T, eq_pred>& up = vector3<T, eq_pred>(0, 1, 0))
		{
			vector3<T, eq_pred> zAxis = norm(at - pos);
			vector3<T, eq_pred> xAxis = norm(cross(up, zAxis));
			vector3<T, eq_pred> yAxis = cross(zAxis, xAxis);

			matrix4x4<T, eq_pred> tmp(xAxis.x,            yAxis.x,            zAxis.x,            0,
									  xAxis.y,            yAxis.y,            zAxis.y,            0,
									  xAxis.z,            yAxis.z,            zAxis.z,            0,
									  -dot(xAxis, pos), -dot(yAxis, pos), -dot(zAxis, pos), 1);

			return tmp;
		}

		/**
		 * Create a right-handed, look-at matrix.
		 *
		 * @param pos   the position of the camera in world space
		 * @param at    the target of the camera in world space
		 * @param up    the up vector of the world, usually (0, 1, 0)
		 *
		 * @return   a right-handed, look-at matrix
		 */
		static matrix4x4<T, eq_pred> lookAtRH(const vector3<T, eq_pred>& pos, const vector3<T, eq_pred>& at, const vector3<T, eq_pred>& up = vector3<T, eq_pred>(0, 1, 0))
		{
			vector3<T, eq_pred> zAxis = norm(pos - at);
			vector3<T, eq_pred> xAxis = norm(cross(up, zAxis));
			vector3<T, eq_pred> yAxis = cross(zAxis, xAxis);

			matrix4x4<T, eq_pred> tmp(xAxis.x,          yAxis.x,          zAxis.x,          0,
									  xAxis.y,          yAxis.y,          zAxis.y,          0,
									  xAxis.z,          yAxis.z,          zAxis.z,          0,
									  -dot(xAxis, pos), -dot(yAxis, pos), -dot(zAxis, pos), 1);

			return tmp;
		}
	};
	///////////////////////////////////////////////////////////////////////////////////////////////




	/**
	 * Calculate the inverse of a 4x4 matrix.
	 */
	template <typename T, typename eq_pred>
	matrix4x4<T, eq_pred> inverse(const matrix4x4<T, eq_pred>& matrix)
	{
		matrix4x4<T, eq_pred> tmp = matrix;
		tmp.inverse();
		return tmp;
	}
	///////////////////////////////////////////////////////////////////////////////////////////////




	/**
	 * Calculate the transpose of a 4x4 matrix.
	 */
	template <typename T, typename eq_pred>
	matrix4x4<T, eq_pred> transpose(const matrix4x4<T, eq_pred>& matrix)
	{
		matrix4x4<T, eq_pred> tmp = matrix;
		tmp.transpose();
		return tmp;
	}
	///////////////////////////////////////////////////////////////////////////////////////////////




	/**
	 * Multiply a 4-dimensional vector with a 4x4 matrix.
	 * The vector is treated as a column and is multiplied by each row of the matrix.
	 *
	 * @return   a 4-dimensional vector that has been transformed by the matrix
	 */
	template <typename T, typename eq_pred>
	vector4<T, eq_pred> operator*(const vector4<T, eq_pred>& l, const matrix4x4<T, eq_pred>& r)
	{
		//         Vector                   Matrix        
		//                          ( 0,0  0,1  0,2  0,3 )
		// ( l.x  l.y  l.z  l.w ) * ( 1,0  1,1  1,2  1,3 )
		//                          ( 2,0  2,1  2,2  2,3 )
		//                          ( 3,0  3,1  3,2  3,3 )
		vector4<T, eq_pred> tmp(dot(l, r.col1()),
								dot(l, r.col2()),
								dot(l, r.col3()),
								dot(l, r.col4()));
		return tmp;
	}
	///////////////////////////////////////////////////////////////////////////////////////////////




	template <typename T, typename eq_pred>
	vector3<T, eq_pred> operator*(const math::vector3<T, eq_pred>& l, const math::matrix4x4<T, eq_pred>& r)
	{
		vector3<T, eq_pred> tmp(l * r.rotation());
		tmp += r.col4().xyz();
		return tmp;
	}
	///////////////////////////////////////////////////////////////////////////////////////////////




	/**
	 * Calculate a hash for a 4x4 matrix.
	 */
	template <typename T, typename eq_pred>
	std::size_t hash_value(const matrix4x4<T, eq_pred>& that)
	{
		std::size_t seed = 0;
		boost::hash_combine(seed, that._11);
		boost::hash_combine(seed, that._12);
		boost::hash_combine(seed, that._13);
		boost::hash_combine(seed, that._14);
		boost::hash_combine(seed, that._21);
		boost::hash_combine(seed, that._22);
		boost::hash_combine(seed, that._23);
		boost::hash_combine(seed, that._24);
		boost::hash_combine(seed, that._31);
		boost::hash_combine(seed, that._32);
		boost::hash_combine(seed, that._33);
		boost::hash_combine(seed, that._34);
		boost::hash_combine(seed, that._41);
		boost::hash_combine(seed, that._42);
		boost::hash_combine(seed, that._43);
		boost::hash_combine(seed, that._44);
		return seed;
	}
	///////////////////////////////////////////////////////////////////////////////////////////////




	/**
	 * In order to compile a binary that actually contains the code for our matrices,
	 * we have to explicitly instantiate the template matrices and export those:
	 */
	template class MATH_API matrix4x4<float>;
	template class MATH_API matrix4x4<double>;

	/**
	 * Because matrix4x4<float> is a real pain to write, we create a typedef.
	 */
	typedef matrix4x4<float>    float4x4;
	typedef matrix4x4<double>   double4x4;
}
///////////////////////////////////////////////////////////////////////////////////////////////////




/**
 * Read a 4x4 matrix from a binary stream.
 */
template <typename T, typename eq_pred>
boost_ext::iostreams::bin_istream& operator>>(boost_ext::iostreams::bin_istream& stream, math::matrix4x4<T, eq_pred>& matrix)
{
	stream >> matrix._11 >> matrix._12 >> matrix._13 >> matrix._14;
	stream >> matrix._21 >> matrix._22 >> matrix._23 >> matrix._24;
	stream >> matrix._31 >> matrix._32 >> matrix._33 >> matrix._34;
	stream >> matrix._41 >> matrix._42 >> matrix._43 >> matrix._44;

	return stream;
}
///////////////////////////////////////////////////////////////////////////////////////////////////

/**
 * Write a 4x4 matrix to a binary stream.
 */
template <typename T, typename eq_pred>
boost_ext::iostreams::bin_ostream& operator<<(boost_ext::iostreams::bin_ostream& stream, const math::matrix4x4<T, eq_pred>& matrix)
{
	stream << matrix._11 << matrix._12 << matrix._13 << matrix._14;
	stream << matrix._21 << matrix._22 << matrix._23 << matrix._24;
	stream << matrix._31 << matrix._32 << matrix._33 << matrix._34;
	stream << matrix._41 << matrix._42 << matrix._43 << matrix._44;

	return stream;
}
///////////////////////////////////////////////////////////////////////////////////////////////////

#pragma pack(pop)
